#########################################################################
# 	SCRIPT FOR "RECONSIDERING THE NEIGHBORHOOD EFFECT"			#
#	Authors: Martin Bisgaard, Peter T. Dinesen & Kim S�nderskov 	#
#	Last edited: 10-11-2015


rm(list=ls())
require(foreign)

#load data
dat<-read.dta("mydirectory")
names(dat)
dim(dat)


####################################
# RESCALING AND RECODING VARIABLES #

require(car)

#main DV: economic perceptions
dat$stfeco<-dat$stfeco/10 #rescaling to 0-1
table(dat$stfeco) 

#self-reported exposure to TV, Radio and newspapers
dat$tvtot<-recode(dat$tvtot,'88=NA;99=NA');dat$tvtot<-dat$tvtot/7
table(dat$tvtot)
dat$rdtot<-recode(dat$rdtot,'88=NA;99=NA');dat$rdtot<-dat$rdtot/7
table(dat$rdtot)
dat$nwsptot<-recode(dat$nwsptot,'88=NA;99=NA');dat$nwsptot<-dat$nwsptot/7
table(dat$nwsptot)

#satisfaction with government (used for supplemental analyses)
dat$stfgov<-dat$stfgov/10
table(dat$stfgov)

#vote intention (used for supplemental analyses)
dat$pvote<-NA
dat$pvote[dat$essround == 1]<-dat$prtvtdk2002[dat$essround == 1]
dat$pvote[dat$essround == 2]<-dat$prtvtdk2004[dat$essround == 2]
dat$pvote[dat$essround == 3]<-dat$prtvtdk2006[dat$essround == 3]
dat$pvote[dat$essround == 4]<-dat$prtvtdk2008[dat$essround == 4]
dat$pvote[dat$essround == 5]<-recode(dat$prtvtdk2010[dat$essround == 5],'4=5;5=6;6=7;7=8;8=12;9=10;10=11')
dat$pvote<-recode(dat$pvote,'66=13;77=13;88=13;99=13')
table(dat$pvote)

#recoding to vote intention to binary: incumbent parties/opposition parties
dat$gov<-NA
dat$gov[dat$pvote%in%c(3,8)]<-1
dat$gov[dat$pvote%in%c(1:2,4:7,9:13)]<-0
table(dat$gov)



# rescaling additional measures
dat$munemploy<-dat$munemploy/100 #municipality unemployment
dat$pincome2000_1k<-dat$pincome2000/100000 
dat$dayslived_y<-dat$dayslived/365 #scaling to years lived at location
dat$lrscale<-dat$lrscale/10 #left-right orientation



################################
#  Code producing Figure 1	#

#Results plotted in Panel A
#computing differences from 80m context

ueb.names<-paste("ueb",c(80,130,180,250,350,500,750,1000,1250,1500,1750,2000,2250,2500),sep="")

diff<-dat[dat$popdens80>9,ueb.names]#matrix with all contextual measures
diff.out<-matrix(NA,nr=nrow(dat[dat$essround!=9&dat$popdens80>9,]),nc=length(ueb.names)-1)
for(i in 1:ncol(diff.out)) diff.out[,i]<-diff[,i+1]-diff[,1]


#Results plotted in Panel B


#inter-correlation between 80m and all other contexts (solid line in figure)

coef.out<-matrix(NA,nr=14,nc=2)#matrix to hold output

for (i in 2:14){
form<-substitute(ueb80 ~ sub,list(sub=as.name(ueb.names[i]))) #model formula
fit1<-lm(form,data=dat[dat$popdens80>9,]) #model fit (restricted to respondents with more than 9 ppl in their 80m radius)
coef.out[i,1]<-coef(fit1)[ueb.names[i]]#extract coefficient
	SSreg<-sum(fit1$residuals^2)/fit1$df
	SStot<-var(fit1$model[,1])
coef.out[i,2]<-1-SSreg/SStot #extract adj. rsquared
}
coef.out
coef.out[1,]<-1 #setting an artificial '1' for 80m on 80m regression
 
#inter-correlation between municipality unemployment and all other contexts (dashed line in figure)
coef.out1<-matrix(NA,nr=14,nc=2)
for (i in 1:14){
form<-substitute(munemploy ~ sub,list(sub=as.name(ueb.names[i])))
fit1<-lm(form,data=dat[dat$essround!=9&dat$popdens80>9,])

coef.out1[i,1]<-coef(fit1)[ueb.names[i]] #coef
	SSreg<-sum(fit1$residuals^2)/fit1$df
	SStot<-var(fit1$model[,1])
coef.out1[i,2]<-1-SSreg/SStot #adj. rsquared
}
coef.out1



#	PLOTTING FIGURE 1	#

#calculating the proportions of distribution below/over 5% (depicted in Panel A, Figure 1)
below=prop.table(table(diff.out[,13]<c(-.05)))[2]
above=prop.table(table(diff.out[,13]>c(.05)))[2]
remain=1-(below+above)
below;above;remain
below+above+remain


#pdf("bivar.pdf",width=7.75,height=4)
setEPS()
postscript("bivar.eps",width=7.75,height=4)

xax<-c(80,130,180,250,350,500,750,1000,1250,1500,1750,2000,2250,2500)#spacing for x-axis (Panel B)

par(mfrow=c(1,2))
dens=density(diff.out[,13],na.rm=T,bw=0.005)
plot(dens,xlim=c(-.3,.2),main="",xlab="",yaxt='n',ylab="Kernel density")
axis(2,las=2)
mtext(1,text="Deviation in unemployment rates",line=2.5)
mtext(1,text="(2,500m - 80m)",line=3.5)
mtext(3,text="A",font=2,cex=2,line=2)

mtext(3,text="How much does aggregate unemployment (2,500m)",line=1,cex=.8)
mtext(3,text="deviate from residential unemployment (80m)?",line=.25,cex=.8)
mtext(2,text="(bandwidth = 0.005)",cex=.8,line=2)


x1<-min(which(dens$x>=.05))
x2<-max(which(dens$x<c(-.05)))

with(dens,polygon(x=c(x[c(x1,x1:x2,x2)]),y=c(0,y[x1:x2],0),col="grey80"))
text(y=0.35,x=-.1,label="19 %")
text(y=0.35,x=0,label="44 %")
text(y=0.35,x=.1,label="37 %")


plot(x=xax,coef.out[,2],ylim=c(0,1),yaxt="n",xaxt="n",pch=1,col="white",
	xlab="",ylab="R-Squared")
mtext(3,text="B",font=2,cex=2,line=2)
mtext(3,text="The proportion of variance explained",line=1,cex=.8)
mtext(3,text="by different contextual measures",line=.25,cex=.8)

lines(xax,coef.out[,2],lwd=1)
lines(xax,coef.out1[,2],lwd=1,lty=2)
axis(2,seq(0,1,.2),las=2,cex.axis=1)

pos<-c(0,250,500,750,1000,1250,1500,1750,2000,2250,2500)
axis(1,at=pos,labels=NA)
text(pos[-1],par()$usr[3]-.06*(par()$usr[4]-par()$usr[3]),labels=c("250m","500m","750m","1000m","1250m","1500m","1750m","2000m","2250m","2500m"),
	adj=1,srt=45,xpd=T,cex=.8)

legend("topright",title="Dependent variable:",bty="n",legend=c("80m unemployment","Municipality unemployment"),lty=c(1,2),cex=.8,xpd=TRUE,lwd=1,pch=c(21,NA,1,NA),pt.bg="black")

points(xax,coef.out[,2],pch=21,cex=1,bg="black")
points(xax,coef.out1[,2],pch=1,cex=1,bg="black")
mtext(1,text="Level of aggregation of contextual",line=2.5)
mtext(1,text="predictor (radius in meters)",line=3.5)

dev.off()




##########################
# Results for Table 1	#
require(memisc)

fit1<-lm(stfeco~factor(municip)+factor(essround)+popdens80+unemplo4+ageR+eduresu+maleR+pincome2000_1k+dayslived_y+lrscale+nwsptot+polintr+ueb80
	,data=dat[dat$popdens80>9,])

fit2<-lm(stfeco~factor(municip)+factor(essround)+popdens80+unemplo4+ageR+eduresu+maleR+pincome2000_1k+dayslived_y+lrscale+nwsptot+polintr+ueb80+ueb2500
	,data=dat[dat$popdens80>9,])

fit3<-lm(stfeco~factor(municip)+factor(essround)+popdens80+unemplo4+ageR+eduresu+maleR+pincome2000_1k+dayslived_y+lrscale+nwsptot+polintr+ueb80+ueb2500+oneunemplo4
	,data=dat[dat$popdens80>9,])

tab<-mtable(fit1,fit2,fit3)

tab<-relabel(tab,
	"ueb80"="Local Unemployment Rate (80m)",
	ueb2500="Local Unemployment Rate (2500m)",
	ueb2500="Local Unemployment Rate (2500m)",
	oneunemplo4="Household Unemployment",
	lrscale="Right--wing",
	polintr="Political interest",
	nwsptot="Newspaper exposure",
	dayslived_y="Years lived at location",
	pincome2000_1k="Income (in 100k DKR)",
	ageR="Age",
	unemplo4="Unemployed",
	popdens80="Population (80m)",
	eduresu="Years of education",
	"maleR: Male/Female"="Male"
	)

tab








######################################
# Results and code for Figure 2	#
# Sensitivity to size of context	#


ueb.names<-paste("ueb",c(80,130,180,250,350,500,750,1000,1250,1500,1750,2000,2250,2500),sep="")
popdens.names<-paste("popdens",c(80,130,180,250,350,500,750,1000,1250,1500,1750,2000,2250,2500),sep="")


eff.out2<-matrix(NA,nr=14,nc=2)
eff.out1<-matrix(NA,nr=14,nc=2)
for (i in 1:14){

lmfit2<-lm(substitute(stfeco~sub1+sub2+factor(municip)+factor(essround)+unemplo4+ageR+eduresu+maleR+pincome2000_1k+dayslived_y+lrscale+nwsptot+polintr+partunemplo4
	,list(sub1=as.name(ueb.names[i]),sub2=as.name(popdens.names[i])))
	,data=dat[dat$popdens80>9,])

eff.out2[i,1]<-coef(lmfit2)[ueb.names[i]]
eff.out2[i,2]<-sqrt(diag(vcov(lmfit2)))[ueb.names[i]]

lmfit1<-lm(substitute(stfeco~sub1+sub2+factor(essround)+unemplo4+ageR+eduresu+maleR+pincome2000_1k+dayslived_y+lrscale+nwsptot+polintr+partunemplo4
	,list(sub1=as.name(ueb.names[i]),sub2=as.name(popdens.names[i])))
	,data=dat[dat$popdens80>9,])

eff.out1[i,1]<-coef(lmfit1)[ueb.names[i]]
eff.out1[i,2]<-sqrt(diag(vcov(lmfit1)))[ueb.names[i]]

}
eff.out2
eff.out1

#weighting estimates with 1st to 99th percentile variation.
w<-matrix(NA,nr=length(ueb.names),nc=2)
for (i in 1:length(ueb.names))w[i,]<-quantile(dat[dat$popdens80>9,ueb.names[i]],c(.01,.99),na.rm=T)
w.out<-apply(w,1,diff)

(eff.out2<-eff.out2*w.out)
(eff.out1<-eff.out1*w.out)



######################################
#	PLOTTING				#
#	With/without MUNICI FEs		#

#pdf('ueb-agg-weighted.pdf',width=6.5,height=4)
setEPS()
postscript("ueb-agg-weighted.eps",width=6.5,height=4)

jit<-c(rep(10,5),15,rep(20,8))
xax<-c(80,130,180,250,350,500,750,1000,1250,1500,1750,2000,2250,2500)

par(mfrow=c(1,1),omd=c(0.05,1,0,1))
plot(xax,eff.out2[,1],main="",yaxt="n",col="white",
	ylim=c(-.1,.05),xaxt="n",
	ylab="Marginal Effect of Local Unemployment Rate",
	xlab="Level of aggregation of contextual predictor (radius in meters)",
	cex.lab=1)
segments(xax-jit,eff.out2[,1],xax-jit,eff.out2[,1]+1.97*eff.out2[,2],lwd=1)
segments(xax-jit,eff.out2[,1],xax-jit,eff.out2[,1]-1.97*eff.out2[,2],lwd=1)
abline(h=0,lty=2)
points(xax-jit,eff.out2[,1],lwd=1,pch=21,bg="black",col="black",cex=1.25)
axis(2,las=2,at=round(seq(-.1,.05,.05),2),cex.axis=1)
legend("bottomleft",lwd=1.5,lty=1,bty="n",legend=c("95% CI"),cex=1)
legend("bottomright",pch=21,pt.bg=c("black","white"),legend=c("With municipality fixed effects","Without municipality fixed effects"),bty="n",
	cex=1,pt.cex=1,col="black")

#ADDING W/O MUNICI FEs
segments(xax+jit,eff.out1[,1],xax+jit,eff.out1[,1]+1.97*eff.out1[,2],lwd=1,col="black")
segments(xax+jit,eff.out1[,1],xax+jit,eff.out1[,1]-1.97*eff.out1[,2],lwd=1,col="black")
abline(h=0,lty=2)
points(xax+jit,eff.out1[,1],lwd=1,pch=21,bg="white",col="black",cex=1.25)


pos<-c(0,250,500,750,1000,1250,1500,1750,2000,2250,2500)
axis(1,at=pos,labels=NA)
text(pos[-1],par()$usr[3]-.06*(par()$usr[4]-par()$usr[3]),labels=c("250m","500m","750m","1000m","1250m","1500m","1750m","2000m","2250m","2500m"),
	adj=1,srt=45,xpd=T,cex=1)

dev.off()





##################################################
# Results and code for Table 2 and Figure 3	#

fit2FE<-lm(stfeco~factor(municip)+factor(essround)+unemplo4+ageR+eduresu+maleR+lrscale+polintr+oneunemplo4
	+pincome2000_1k+dayslived_y+popdens80+ueb80*nwsptot+munemploy*nwsptot
	,data=dat[dat$popdens80>9,])

fit1FE<-lm(stfeco~factor(municip)+factor(essround)+unemplo4+ageR+eduresu+maleR+lrscale+polintr+oneunemplo4
	+pincome2000_1k+dayslived_y+popdens80+ueb80*nwsptot+popdens2500+ueb2500*nwsptot
	,data=dat[rownames(dat)%in%rownames(fit2FE$model),])

tab<-mtable(fit1FE,fit2FE)

tab<-relabel(tab,
	"ueb80"="Local Unemployment Rate (80m)",
	ueb2500="Local Unemployment Rate (2500m)",
	ueb2500="Local Unemployment Rate (2500m)",
	oneunemplo4="Household Unemployment",
	lrscale="Right--wing",
	polintr="Political interest",
	nwsptot="Newspaper exposure",
	dayslived_y="Years lived at location",
	pincome2000_1k="Income (in 100k DKR)",
	ageR="Age",
	unemplo4="Unemployed",
	popdens80="Population (80m)",
	eduresu="Years of education",
	"maleR: Male/Female"="Male"
	)

tab

#F-tests for interaction terms in FE models
fit20<-lm(stfeco~factor(municip)+factor(essround)+unemplo5+ageR+eduresu+maleR+lrscale+polintr
	+pincome2000_1k+dayslived_y+popdens80+ueb80+munemploy+nwsptot
	,data=dat[rownames(dat)%in%rownames(fit2FE$model),])
fit0<-lm(stfeco~factor(municip)+factor(essround)+unemplo5+ageR+eduresu+maleR+lrscale+polintr
	+pincome2000_1k+dayslived_y+popdens80+ueb80+nwsptot+popdens2500+ueb2500
	,data=dat[rownames(dat)%in%rownames(fit2FE$model),])
anova(fit0,fit1FE)
anova(fit20,fit2FE)


######################################
# Code for Figure 3			#

summary(dat$munemploy)
news<-seq(0,1,length.out=10)
muni<-seq(min(dat$munemploy,na.rm=TRUE),max(dat$munemploy,na.rm=TRUE),length.out=10)
ueb25<-seq(min(dat$ueb2500,na.rm=TRUE),max(dat$ueb2500,na.rm=TRUE),length.out=10)

fit1<-fit1FE #figure was originally without municipality FEs and this adjusts the code to models /w FEs
fit2<-fit2FE


#coefs
b1<-fit1$coef["ueb2500"]
b3<-fit1$coef["nwsptot:ueb2500"]
#SEs
v1<-diag(vcov(fit1))["ueb2500"]
v3<-diag(vcov(fit1))["nwsptot:ueb2500"]
cov3<-vcov(fit1)["nwsptot:ueb2500","ueb2500"]

se<-rep(NA,10)
cf<-rep(NA,10)
for(i in 1:10){
	se[i]<-sqrt(v1+news[i]^2*v3+2*news[i]*cov3)
	cf[i]<-b1+b3*news[i]
	}
se
cf

w<-matrix(NA,nr=length(ueb.names),nc=2)#weighting wrt to variation (1st to 99th percentile)
for (i in 1:length(ueb.names))w[i,]<-quantile(dat[dat$essround!=9&dat$popdens80>9,ueb.names[i]],c(.01,.99),na.rm=T)
w.out<-apply(w,1,diff)

se<-se*w.out[14]	
cf<-cf*w.out[14]


#coefs
b1<-fit2$coef["munemploy"]
b3<-fit2$coef["nwsptot:munemploy"]
#SEs
v1<-diag(vcov(fit2))["munemploy"]
v3<-diag(vcov(fit2))["nwsptot:munemploy"]
cov3<-vcov(fit2)["nwsptot:munemploy","munemploy"]

se1<-rep(NA,10)
cf1<-rep(NA,10)
for(i in 1:10){
	se1[i]<-sqrt(v1+news[i]^2*v3+2*news[i]*cov3)
	cf1[i]<-b1+b3*news[i]
	}
se1
cf1

w.mun<-quantile(dat[dat$essround!=9&dat$popdens80>9,"munemploy"],c(0.01,.99),na.rm=T)
w.mun<-w.mun[2]-w.mun[1]

se1<-se1*w.mun
cf1<-cf1*w.mun

#coefs
b1<-fit2$coef["ueb80"]
b3<-fit2$coef["ueb80:nwsptot"]
#SEs
v1<-diag(vcov(fit2))["ueb80"]
v3<-diag(vcov(fit2))["ueb80:nwsptot"]
cov3<-vcov(fit2)["ueb80:nwsptot","ueb80"]

se2<-rep(NA,10)
cf2<-rep(NA,10)
for(i in 1:10){
	se2[i]<-sqrt(v1+news[i]^2*v3+2*news[i]*cov3)
	cf2[i]<-b1+b3*news[i]
	}
se2
cf2


se2<-w.out[1]*se2
cf2<-w.out[1]*cf2




#######PLOTTING####
#pdf('nwsp-int.pdf',width=8,height=3)
setEPS()
postscript('nwsp-int.eps',width=8,height=3)


par(mfrow=c(1,3),omd=c(.05,1,0,1))
#par(mfrow=c(1,1))
plot(news,cf2,type="l",ylim=c(-.2,.1),cex.axis=.8,,main="",
	ylab="",xaxt="n",yaxt="n",
	xlab="")
axis(2,at=seq(-.2,.1,.05),las=2,cex.axis=1)
axis(1,cex.axis=1)
mtext(2,text="Marginal effect of unemployment rate",line=5)
mtext(2,text="(1st to 99th percentile)",line=3.75,cex=.8)
mtext(3,text="A",font=2,line=2,cex=1.5)
mtext(3,text="Unemployment in 80m radius",font=1,line=0.5,cex=.8)



abline(h=0,lty=2)
#X.vec<-c(news,1,rev(news),0)
#Y.vec<-c(cf-1.97*se,cf[10]+1.97*se[10],rev(cf+1.97*se),cf[1]+1.97*se[1])
#polygon(X.vec,Y.vec,col="grey90",border=NA);lines(news,cf,col="black")

lines(news,cf2+1.97*se2)
lines(news,cf2-1.97*se2)

par(new=TRUE)
plot.new()
hist(dat$nwsptot[dat$popdens80>9&dat$essround!=9],xpd=TRUE,col="grey95",breaks=na.omit(unique(dat$nwsptot)),ylim=c(0,15),
	freq=FALSE,axes=FALSE,xlab="",ylab="",main="")

#dev.new(width=4,height=4.7)
#par(mfrow=c(1,1))
plot(news,cf,type="l",ylim=c(-.2,.1),cex.axis=.8,,main="",
	ylab="",xaxt="n",yaxt="n",
	xlab="")
axis(2,at=seq(-.2,.1,.05),las=2,cex.axis=1)
axis(1,cex.axis=1)

mtext(3,text="B",font=2,line=2,cex=1.5)
mtext(3,text="Unemployment in 2500m radius",font=1,line=0.5,cex=.8)
mtext(1,text="Newspaper exposure",line=3)
abline(h=0,lty=2)
#X.vec<-c(news,1,rev(news),0)
#Y.vec<-c(cf-1.97*se,cf[10]+1.97*se[10],rev(cf+1.97*se),cf[1]+1.97*se[1])
#polygon(X.vec,Y.vec,col="grey90",border=NA);lines(news,cf,col="black")

lines(news,cf+1.97*se)
lines(news,cf-1.97*se)

par(new=TRUE)
plot.new()
hist(dat$nwsptot[dat$popdens80>9&dat$essround!=9],xpd=TRUE,col="grey95",breaks=na.omit(unique(dat$nwsptot)),ylim=c(0,15),
	freq=FALSE,axes=FALSE,xlab="",ylab="",main="")


#dev.new(width=4,height=4.7)
#par(fig=c(0,1,0,1))
plot(news,cf1,type="l",ylim=c(-.2,.1),cex.axis=.8,main="",
	ylab="",xaxt="n",yaxt="n",
	xlab="")
axis(2,at=seq(-.2,.1,.05),las=2,cex.axis=1)
axis(1,cex.axis=1)

mtext(3,text="C",font=2,line=2,cex=1.5)
mtext(3,text="Unemployment in municipality",font=1,line=0.5,cex=.8)

lines(news,cf1+1.97*se1)
lines(news,cf1-1.97*se1)
abline(h=0,lty=2)

par(new=TRUE)
plot.new()
hist(dat$nwsptot[dat$popdens80>9&dat$essround!=9],xpd=TRUE,col="grey95",breaks=na.omit(unique(dat$nwsptot)),ylim=c(0,15),
	freq=FALSE,axes=FALSE,xlab="",ylab="",main="")

dev.off()




#END OF SCRIPT
